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A connection is established between discrete stochastic model describing microscopic motion of 
fluctuating cells, and macroscopic equations describing dynamics of cellular density. Cells move to- 
wards chemical gradient (process called chemotaxis) with their shapes randomly fluctuating. Non- 
linear diffusion equation is derived from microscopic dynamics in dimensions one and two using 
excluded volume approach. Nonlinear diffusion coefficient depends on cellular volume fraction and 
it is demonstrated to prevent collapse of cellular density. A very good agreement is shown be- 
tween Monte Carlo simulations of the microscopic Cellular Potts Model and numerical solutions 
of the macroscopic equations for relatively large cellular volume fractions. Combination of micro- 
scopic and macroscopic models were used to simulate growth of structures similar to early vascular 
networks. 
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I. INTRODUCTION 

So far most models used in biology have been de- 
veloped at specific scales. Establishing a connection 
between discrete stochastic microscopic description and 
continuous deterministic macroscopic description of the 
same biological phenomenon would allow one to switch 
when needed from one scale to another, considering 
events at individual (microscopic) cell level such as cell- 
cell interaction or cell division to events involving thou- 
sands of cells such as organ formation and development. 
Due to the fast calculation speed possible with the contin- 
uous model, one can quickly test wide parameter ranges 
and determine satiability conditions and then use this 
information for running Monte Carlo simulations of the 
stochastic discrete dynamical systems. Also, continuous 
models provide very good approximation for systems con- 
taining a biologically reahstic (i.e., large) number of cells, 
for which numerical simulations of stochastic trajectories 
can be prohibitive. 

Most continuous biological models have been postu- 
lated either by requiring certain biologically relevant fea- 
tures from the solutions or making it easier to analyze 
behavior of solutions using certain mathematical tech- 
niques. In particular, system of nonlinear partial differ- 
ential equations model with chemotactic term was used 
in P, Q to simulate the de novo blood vessel formation 
from the mesoderm. The rational for the model was pro- 
vided by the experimental observations 0] demonstrating 
that chemotaxis played an important role in guiding cells 
during early vascular network formation. Discrete mod- 
els have been also applied to simulating vasculogenesis 
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and angeogenesis 0, H, Q . 

In this paper we derive continuous macroscopic limits 
of the ID and 2D microscopic cell-based models with ex- 
tended cell representations, in the form of nonlinear dif- 
fusion equations coupled with chemotaxis equation. We 
demonstrate that combination of the discrete model and 
derived continuous model can be used for simulating bio- 
logical phenomena in which a nonconfluent population of 
cells interact directly and via diffusible factors, forming 
an open network structures in a way similar to forma- 
tion of networks during vasculogenesis P, i^] and pattern 
formation in limb cell cultures 7]. 

Continuous limits of microscopic models of biological 
systems based on point-wise based cell representation 
were extensively studied over the last 30 years. The 
classical Keller-Segel PDE model has been derived in 
Q from a discrete model with point-wise cells under- 
going random walk in chemotactic field and then studied 
in [13, [O, [13] ■ Cells in this model secret diffusing 
chemical at constant rate and detect local concentration 
c of this chemical due to process called chemotaxis. The 
chemical is called an attractant or repellent depending 
whether cell moves towards chemical gradient or in op- 
posite direction. Aggregation occurs if attraction exceeds 
diffusion of cells. For point-wise cells aggregation results 
in infinite cellular density corresponding to the solution 
of the macroscopic Keller-Segel equation becoming infi- 
nite in finite time (also called blow up in finite time or 
collapse of the solution) [H, • 

There have been few attempts to derive macroscopic 
limits of microscopic models which treat cells as extended 
objects consisting of several points. In [l^ the diffusion 
coefficient for a collection of noninteracting randomly 
moving cells was derived from a one-dimensional Cel- 
lular Potts Model (CPM). A microscopic limit of sub- 
cellular elements model [16] was derived in the form of 
an advection-diffusion PDE for cellular density. In our 
previous papers [13, [l8|, [l9| we studied the continuous 
limit of the CPM describing individual cell motion in 



a medium and in the presence of an external field witl 
contact cell-cell interactions in mean-field approximation 
However, mean-field approximation does not allow one tc 
consider high density of cells when cellular volume frac 
tion (fraction of volume occupied by cells) (/) is of th( 
order one. 

In this paper we go beyond mean-field approximation 
Namely, we take into account finite size of cells in th( 
CPM, use exclusion volume principle (meaning that tw( 
cells can not occupy the same volume) and derive the fol 
lowing macroscopic nonlinear diffusion equation for evo 
lution of cellular density p{r,t) in one (ID) 
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dtP - 

and two dimensions (2D): 
dtp 



1 



XoVr- [pVrc(r,i)], (1 



.l-0 + (/.log(0) 

-XoVr- [pVrc(r,i)], {2 

which do not have blow up in finite time. These nonlin 
ear diffusion equations are coupled with the equation fo: 
evolution of chemical field c(r, t) : 

dtc{r,t) = DcVlc + ap-"fc. (3 

Here (/) is a volume fraction (fraction of volume occupiec 
by cells) . In ID case cells have a form of fiuctuating rod: 

and 4> = L^x^p(r,t), where L^'' is an average length Oi 
cells. In 2D case we assume that cells are fluctuating 



rectangles and (f> = L^x^ L^^ p{v,t)^ where (Vx' , Ly'') are 
average length and width of a cell. Here p{r, t) is the den- 
sity of cells normalized by the total number of cells N : 
J p{r, t)dr = N and r is a vector of spatial coordinates in 
ID or 2D. D2 is the diffusion coefficient for a motion of 
an isolated cell, xo defines strength of chemotactic inter- 
actions, Dc is the diffusion coefficient of a chemical, a is 
the production rate of a chemical and 7 is the decay rate 
of a chemical. Typical microscopic picture of distribu- 
tion of individual cells is shown in Figure [1] for 2D case. 
Solutions of Eqs. (H]), ^ and ^ describe coarse-grained 
macroscopic cellular dynamics. Below we show very good 
agreement even for relatively large densities 4> ~ 0.3. , be- 
tween solutions of Eqs. ((T|) , ([2]) , ((31) and ensemble average 
of stochastic trajectories of the microscopic CPM. 

In Section II we introduce general microscopic equa- 
tions describing motion of cells based on random fluc- 
tuations of their shapes and their interactions. We as- 
sume that each cell has a rectangular shape and consider 
stochastic differential equations for motion of cells as well 
as Smoluchowski equation (see e.g. 2Ql) for multi-cellular 
probability density function. In Section HI we consider 
a particular case of microscopic cellular dynamics rep- 
resented by the CPM without excluded volume interac- 
tions. Coefficients of Smoluchowski Eq. are derived from 
the CPM and stochastic dynamics of shapes and posi- 
tions of cells are reduced to solutions of the closed equa- 
tion describing positions of cells. In Section IV we con- 
sider cell motion and cell-cell interactions with collisions 
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FIG. 1: Representation of cells in the 2D CPM in the form of 
fluctuating rectangles. 539 cells (15% volume fraction) with 



r(0) 



-(0)x 



1.666667 are shown. 



resolved through the jump processes resulting in equa- 
tions ^ , ^ and ^ . In Section V numerics for the con- 
tinuous macroscopic equations is compared with Monte 
Carlo simulations of the microscopic CPM. In Section VI 
we summarize main results and discuss future directions. 



II. MICROSCOPIC MOTION OF CELLS 

Motion of many eukaryotic cells and bacteria is ac- 
companied by the random fluctuations of their shapes 
pll . 12^ . l23j resulting in a random diffusion of center of 
mass of an isolated cell. Coefficient of such diffusion can 
be measured experimentally (see e.g. IM])- Fluctuations 
of cellular membrane in the presence of a chemical field 
are more likely in the direction of chemical gradient (for 
chemoattractant) of in opposite direction (for chemore- 
pellent). Cells can also interact through direct contact 
which includes cell-cell adhesion and can be modeled us- 
ing excluded volume principle. Cellular environment is 
highly viscous and inertia of cell can be ignored. 

In this paper we assume that each cell has a fluctu- 
ating rectangular shape and allow random fluctuation 
of the dimensions of each rectangular cell (see Figures 
[T] and [2]) . Positions and shapes of cells are completely 
characterized by a finite set of dynamical variables in the 
configuration space: X = (Ri, Li, R2, L2, . . . , Rat, Lat), 
where N is the total number of cells, Rj is a position 
of center of mass of jth cell, hj is the size of jth cell in 
D spatial dimensions. We consider D = 1 and D = 2 
in which case cells are moving over substrate but results 
can be extended to the D = 3 case. Microscopic descrip- 
tion is provided by the multi-cellular probability density 
function (PDF) P{X,t) defined as ensemble average (•) 



over stochastic trajectories (R' (t), L' (t)), j = 1, , 



,N 



in the configuration space X: P{X,t) = (O^i '^(■^^ 
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FIG. 2; 2D CPM cell representation. Grey and white colors 
are used to indicate cell body and surrounding extracellular 
matrix respectively. Cell can grow or shrink in x and y di- 
rection by adding or removing one row (or column) of pixels, 
pixels. 



'R!j{t))5{Lj — 'L'j{t))) determined by solution of coupled 
stochastic equation 



^ = A(X,c,t) + B(X,c,<)e(t), 



(4) 



and chemotaxis equation for the chemical field c(r, i). 
Here r is the spatial coordinate, A(X, t) is 2DN- 
component vector, I3{X,t) is 2DN x 2DN matrix and 
^(t) is 2£'A^-component stationary Gaussian stochastic 
process with zero correlation time and zero mean 

(0>=0, m)^jit'))^6,Jit-t'), i,j = l,...,2Z?7V,(5) 

where 6i_j is a Kronecker's symbol. 

Application of the Stratonovich stochastic calculus to 
the Eq. (HI) results in a "multi-cellular" Fokker-Planck 
equation [20] in configuration space X 



dtP{X,t) 



-V • [vP] + V • [L» • VP] , 

V = A--B -W ■ B^, 
2 

1 . 



(6) 



where D{X,t) is a 2DN x 2DN diffusion matrix, V = 
(9ri , 9li , • ■ • , i9rn : ^w) is the gradient operator in 2DN 
dimensions and 13^ represents transposed matrix B. 

Dynamics of the chemical field c(r, t) is described by a 
diffusion equation 

5tc(r, t) = V,D, • Vrc + j A{X, r, t)P{X, t)dX - -fc,{7) 

where Vr is the gradient operator for the spatial co- 
ordinate and diffusion coefficient for the chemical field 
Dc in general case can depend on r and t. Term 



/ A{X^Y^t)P{X,t)dX describes production of chemical 
by cells at the rate of A{X, r, t) and 7 is a decay rate of 
the chemical. Cells produce chemical into intercellular 
space through their membranes. Therefore, in 3D case 
A{X, r, t) is nonzero only at the cellular membrane. How- 
ever, in ID and 2D cells are moving on substrate, while 
chemical diffuse over entire three dimensional space so in 
ID and 2D cases A{X, r, t) is nonzero inside cells. 

We assume that v has a form of a potential v — 
— /3V4>(X, t), where (3 = 1/T is the inverse effective tem- 
perature T of the chemical fluctuations. Multi-cellular 
Fokker-Planck equation ^ is then reduced to multi- 
cellular Smoluchowski equation: 



dtP{X, t) = V -b- VP + /3(V$)P , 



(8) 



(More general case of v being any function can be stud- 
ied using similar approach.) If we neglect fluctuations 
of the cellular size Lj and chemotaxis, c = 0, then Eq. 

is similar to the Smoluchowski Eq. for the Brownian 
dynamics of colloidal particles (see e.g. [2^ for a review) 
and Eq. (j4]) has a form of the Langevin equation for 
interacting Brownian particles with term B£^ represent- 
ing thermal forces from solution in colloids. However, in 
general case considered here both chemotaxis and fluctu- 
ations of cellular shape are taken into account. Mecha- 
nisms of random fluctuation of cellular shape and cellular 
motion are still not completely clear and subject of active 
research [HI, [13, [H . 

We impose excluded volume constraint by choosing 
f^{X,c) = 00 if any two cells overlap. We assume in 
what follows that all direct interactions between cells is 
of this type. We also allow indirect interactions between 
cells mediated by chemotaxis for which we choose 



$(X,c) 
$(X,c) 



$(X,c) 



TV 

- Ex(L,Mr,t) 

i=i 

if cells do not overlap, 

= 00 

if any pair of cells overlaps. 



where x(Lj) represents strength of chemotactic interac- 
tion as a function of cellular sizes Lj and <&(X, c) rep- 

resents chemotaxis-independent terms of the potential. 
We assume for simplicity that chemotactic interaction 
depends only on gradient of c(r, t) at the center of mass 



of each cefl. $(X, ( 



is also responsible for preserving 



cellular shape close to some equilibrium shape. Without 
this term shape (size) of each cell would experience un- 
bounded random fluctuation which is non biological. We 
consider specific form of $(X, c) in Section UlI CI 

Our main goal is to derive a macroscopic equation de- 
scribing dynamics of (total) cell probability density func- 
tion 



P(r, t) 



N 



(9) 



4 



coupled with c(r,t). from microscopic equations ([7]) and 

N N 

m. Herepj(r,i) = J P{X,t) JJ dRi JJ dL„ 

1 = 1, l^j m=l Rj=r 

is a single-cell probability density function of the posi- 
tion of center of mass. After approximating A{X, r, t) = 
a^^-^ <5(r — Rj), a = const and assuming that I?c = 
const, Eq. ([7]) is reduced with the help of Eq. © to 
Eq. ([3]). This approximation is justified because typical 
diffusion of a chemical is much faster than cell diffusion 
D2/DC < 1. E.g. D2/DC ~ 1/40 - 1/400 for the cellu- 
lar slime mold Dictyostelium ]2 &, D 2IDC ^ 1 /30 and for 
microglia cells and neutrophils [271 [28| . 



III. MICROSCOPIC CELLULAR DYNAMICS IN 
CELLULAR POTTS MODEL 

Stochastic discrete models are used in a variety of prob- 
lems dealing with biological complexity. One motivation 
for this approach is the enormous range of length scales 
of typical biological phenomena. Treating cells as simpli- 
fied interacting agents, one can simulate the interactions 
of tens of thousands to millions of cells and still have 
within reach the smaller-scale structures of tissues and 
organs that would be ignored in continuum (e.g., par- 
tial differential equation) approaches. At the same time, 
discrete stochastic models including the Cellular Potts 
Model (CPM) can be made sophisticated enough to re- 
produce almost all commonly observed types of cell be- 
havior [2i,[33,[3l|,[33,|3l,[33. Recent book [H reviews 
many of the cell-based models. 

The cell-based stochastic discrete CPM, which is an 
extension of the Potts model from statistical physics, has 
become a common technique for simulating complex bio- 
logical problems including embryonic vertebrate limb de- 
velopment d^, , tumor growth ^37j] and vasculogenesis 
The CPM can be made sophisticated enough to re- 
produce almost all commonly observed types of cell be- 
havior. It consists of a list of biological cells with each cell 
represented by several pixels, a list of generalized cells, 
a set of chemical diffusants and a description of their bi- 
ological and physical behaviors and interactions embod- 
ied in the effective energy E, with additional terms to 
describe absorption and secretion of diffusants and ex- 
tracellular materials. Distribution of multidimensional 
indices associated with lattice cites determines current 
cell system configuration. The effective energy of the sys- 
tem, E, mixes true energies, like cell-cell adhesion, and 
terms that mimic energies, e.g., volume constraint and 
the response of a cell to a gradient of an external field 
(including chemotactic filed) and area constraint. 



A. Cellular Potts Model for cells of rectangular 
shape 

For simplicity, we use in this paper CPM with rectan- 
gular cellular shapes. We also assume that all cells have 



the same type. The results can be extended to the general 
case of the CPM with arbitrary cellular shapes. Also, the 
approach is not limited to using CPM. For example, one 
could use microscopic off-lattice models [l^. 111], where 
each cell is represented by a collection of subcellular ele- 
ments with postulated interactions between them. 

Notice that reduced representation of cells as fluctuat- 
ing rectangles corresponds to intermediate level of de- 
scription where fluctuations of cellular shapes are re- 
placed by fluctuations of cellular sizes. Stochastic Eq. 
@ or Smoluchowski Eq. ([5]) coupled with ^ can be 
used for modeling cell agregation. 

In the CPM change of a cell shape evolves according 
to the classical Metropolis algorithm based on Boltzmann 
statistics and the following effective energy 



E ^E 



ECM 



Ea 



dhe 



H~ -^Perimeter ^" 



lid- 



(10) 



If an attempt to change index of a pixel in a cell leads 
to energy change Ai?, it is accepted with the probability 



$(A£;) 



1, 

-/3AS 



A£; < 
A£; > 0, 



(11) 



where l/j3 — T represents an effective boundary fluctu- 
ation amplitude of model cells in units of energy. Since 
the cells' environment is highly viscous, cells move to 
minimize their total energy consistent with imposed con- 
straints and boundary conditions. If a change of a ran- 
domly chosen pixels' index causes cell-cell overlap it is 
abandoned. Otherwise, the acceptance probability is cal- 
culated using the corresponding energy change. The ac- 
cepted pixel change attempt results in changing location 
of the center of mass and dimensions of the cell. 

We consider 2D case with rectangular shape of each 
cell with sizes Lj = (L^j, Lyj) and position of cen- 
ter of cellular mass at Rj = {xj,yj). Cell motion 
and changing shape are implemented by adding or re- 
moving a row or column of pixels (see Figure We 
assume that cells can come into direct contact and 
that they interact over long distances through chemo- 
taxis. Term Eecm in the Hamiltonian pop phenomeno- 
logically describes net adhesion or repulsion between 
the cell surface and surrounding extracellular matrix: 

Eecm = SjLi "^JECMiL^.j + Lyj), where Jecm is the 
binding energy per unit length of an interface. Term 
EAdhesion = JaLcontact in the Hamiltoniau (HH]) corre- 
sponds to the cell-cell adhesion, where Ja is the binding 
energy per unit length of an interface and Lcontact is the 
total contact area between cells. Term Ep^rimeter de- 
fines an energy penalty function for dimensions of a cell 
deviating from the target values Lt^^^^ : Eperimeter = 

^2^=1 K{Lx,j - LtJ"^ + \y(Ly^j - Lt^Y where A^; and 
\y are Lagrange multipliers. Cells can move up or down 
gradients of both diffusible chemical signals (chemotaxis) 
and insoluble ECM molecules (haptotaxis) described by 
Epieid = SjLi M c(Rj, f)La; jLyj , whcrc /i is an effective 
chemical potential. 



5 



In this paper we neglect cell-cell adhesion Ja = and 
the Hamiltonian pU|) is reduced in 2D to the following 
expression 



JV 



E{X,t) ^^E{r,L,t) 
i=i 



r— R,, . L— 



E{r, L, t) = 2J„r.{Lx + Ly) + K{Lx - Lr^f 

+Xy{Ly~LT^f +fic{r)LxLy. (12) 



B. Master equation for discrete cellular dynamics 

We now represent CPM dynamics by using P(r,L,t), 
probability density for any cell with its center of mass 
at r to have dimensions L at time t. Which means that 
we consider one-cell PDF P(r, L, t) rather than A^-cell 
PDF P{X,t). Let eAr x eAr be the size of a lattice 
site with e <C 1 and let vectors ei^2 indicate changes in 
X and y dimensions: ei — Ar(l,0), 62 = Ar(0, 1). We 
normalize the total probability to the number of cells: 
/ P{r, L, t)drd'L = N. P(r, L, t) used below should not be 
confused with multi-cellular probability density P{X,t). 
Those depend on different arguments. Excluded vol- 
ume constraint implies that position r' and size L' of 
any neighboring cell should satisfy 2\x — x'\ > + L'^, 
2\y-y'\>Ly + L'y. 

Discrete stochastic cellular dynamics under these con- 
ditions is described by the following master equation [l^ : 



P(r,L,t + At) = EL{[i-%,/(r 



ee^-; r, L, t) 



— r2j,r(r + f L -I- ee^; r, L, t) — i^j,i{r + f e^, L 
— eej; r, L, t) — fij,r(r — f e^, L — ee^; r, L, i)] P(r. L, t) 
+Qj i(r, L; r -f fej, L - eej,t)P(r -I- fe^-, L - ee^, i) . „x 
-H%^(r,L;r- fe,,L-ee,,i)P(r- §e,,L--- +^ ^ 



+fljj{r,L;r - |ej,L + eej,t)P{v ^ 
-H%r(r,L;r + |ej,L + eej,i)P(r-f 



ee 



J' 



2 

-e 

2 J' ■ 



t) 
t) 

eej,t)}. 



ee 



'3 ' 



We incorporate dynamics into MC algorithm by defin- 
ing the time step At. Individual biological cells ex- 
perience diffusion through random fluctuations of their 
shapes. Diffusive coefficient can be measured experimen- 
tally (see e.g. [IS)- We choose At to match experi- 
mental diffusion coefficient. Ea. (|13p would determine a 
version of kinetic/dynamic MC algorithm (see e.g. [s^) 
if Ai were to be allowed to fluctuate. For simplicity we 
assume that At = const. Also i7j.;(r, L; r', L', i) and 
rij_r(r, L; r', L', t) denote probabilities of transitions from 
a cell of length L' and center of mass at r' to a cell 
of dimensions L and center of mass at r. Subscripts 
I and r correspond to transitions by addition/removal 
of a row/colomn of pixels from the rear/lower and 
front/upper ends of a cell respectively. 

We define O^- (r, L; r', L') ee r,(,)(r, L; r', L')[1 - 
0i,i(r)(r, L,t)] where r;(r)(r, L; r', L', t) denote proba- 
bilities of transitions from a cell of length L' and cen- 
ter of mass at r' to a cell of dimensions L and center 



of mass at r without taking into account excluded vol- 
ume principle and cell-cell adhesion. According to the 
CPM we have that r,(r, L; r', L') = r^(r, L; r', L') = 

l<^(^E{r,L) - E{r',V)^ where the factor of 1/8 is due 

to the fact that there are potentially 8 possibilities for 
increasing or decreasing of and Ly (it means that we 
can add/romove pixel from any 4 sides of rectangular 
cell). Second term (f)ji(^^-^{r, L,t) takes into account con- 
tact interactions between cells. It includes contributions 
from 3 possible types of stochastic jump processes due to 
contact interactions between cells: (a) a cell adheres to 
another one, (b) two adhered cells dissociate from each 
other due to membrane fluctuations, (c) membranes of 
two adhered cells are prevented from moving inside each 
other (due to excluded volume constraint) resulting in a 
negative sign of a contribution to a jump probability. If 
neither of these 3 processes happens at given time step 
then 0j ;(r)(r, L,t) = 0. 



Macroscopic limit of master equation and 
mean-field approximation 



Eq. (|T3l) is however not closed because one has 
to know multi-cellular probability density to determine 
<^jj{r){^, L, i). One could use BBGKY-type hierarchy 
[40| similar to the one used in kinetic theory of gases, 
which expresses iteratively n-cell PDF through n + 1-cell 
PDF with truncation at some order. This is however 
extremely difficult and ineffective for large n. Instead, 
we develop in Section IV a nonperturbative approach to 
derivation of Eqs. ([T]) and Q. 

In previous work [TtI . nal . [l^ we studied a macro- 
scopic limit e <C 1 of master Eq. (fT3)) by both ne- 
glecting contact interactions between cells [U [3 ^^'^ 
including contact interactions in mean-field approxima- 
tion [1^, which represents simplest closure for BBGKY- 
type hierarchy. If contact interactions are neglected, then 
(jliLj_t) = and yields in macroscopic limit 

e< 1 u^,m- 

dtP{r, L, t) = D^iVl + AVl)P + 8D2f3XxdL^ {L,P) 
+8D2(3XydL^ [LyP) + D2l3L,Lyfidr [PVrc] , 



A. 



[Jcni + K{Lx - LtJ + -Lyfj.c{r)] , 



Ay I 



Do = 



{Arf 
16At 



Similar Eq. in ID case (with only coordinate x and length 
Lx present) was obtained in Ref. [l3|- 

Mean field approximation assumes decoupling of multi- 
cellular PDF, 



JV 



P{X,t) 



-N 



Y[P{r,L,A) 



r=R, 



(15) 



6 



where factor is due to an assumed normalization 

/ P{r,'L,t)drd'L = N. Mean field approximation is ex- 
act if contact interactions are neglected. This holds since 
we assumed above that chemotaxis depended on aver- 
aged density ([9]) only. Eq. (fTS]) results in decoupling of 
Eq. IH]) into independent Eqs. for P(Rj,Lj,i) for all 
j. This allows direct comparison of Eq. (|14p with Eq. 

and yields that diffusion matrix has a diagonal form 
with main diagonal £>2(1, 1, 4, 4, 1, 1, 4, 4 . . .) in 2D. For 
ID and 3D cases numbers 1 and 4 repeat themselves with 
period 1 and 3, respectively. 

Further comparison of ([5]), ([H]) and (|14p leads to ex- 
pression x(Lj) = Lx.jLyj^ and 



^X,t) ^ E{X,t), 



(16) 



where E{X,t) is given by P^ . 

Taking into account contact cell-cell interaction (ex- 
cluded volume) yields that Eq. ifTSj) is not exact any 
more. Also, potential ^{X,t) is given by (|16p only if 
cells do not overlap. Otherwise (^[X,t) = oo according 

to m. 



Here £'™i„ = £'(r,L(""")) is the minimal value of (fT2l) 



achieved at L = L^""'); 



L 



E 



: £;(r,L("""'), 
XyLT^)^lc{r) 



T {min) 
V 



^2c(r)2 

-^Xx{Jcm - XyLTy) + 2{Jcm - X^Lt^) ^c{r) 



AXxXy 



^^c(r)2 



(20) 



and Z{Y.t) 

2tt 



= (2eAr)2^Lexp(-/3ASze«st/«) ^ 
, is an asymptotic formula for a partition 

function as e ^ 0. (See [l^ for details about convergence 
rate for the case without contact interactions) 

Also, under these conditions we can use Eq. (fT7|) to 
integrate Eq. over L which results in the following 
evolution equation for the the cellular probability density 
p(r,t) = /(P(r,L,i)dL [13, 111: 



dtp = D2VIP - xo Vr • [p Vrc(r, t)] , 



(21) 



D. Boltzmann-like distribution and macroscopic 
equation for cellular density 

From Eqs. (dU, ^ and ^ it follows that typ- 
ical fluctuations of cell sizes are determined by 
PXx(y)L'^i^y-j ^ 1. Suppose xq and t/o are typical scales 
of P(r, L, t) with respect to x and y. We assume that 
(3xqXx » 1 and PuoXy » 1, meaning that xq ^ 
Lx, Uo ^ Ly. We also assume that chemical field c(r,t) 
is a slowly varying function of r on the scale of typical cell 
length meaning that Xc/L^ 3> 1, Vc/Ly 3> 1, where Xc 
and Uc are typical scales for variation of c(r, t) in x and 
y. We also make an additional biologically relevant as- 
sumption that AXxXy ^ M^c(r, i)^ meaning that change 



(chemo) 



is small 



of typical cell size due to chemotaxis 5L 

J {chemo) ^ j (min) 
\^^x(y) I < ^x(y) ■ 

If all these conditions are satisfied then we found by 
using both solutions of Eq. and MC simulations of 
CPM with general initial conditions, that PDF P(r, L, t) 
quickly converges in t at each spatial point r to the fol- 
lowing Boltzmann - like form 



P(r,L,i) =PBo;t^(r,L)p(r,t), 



(17) 



where 



PBo;tz(r,L) 



Z{v) 



eM-P^Ei^ngth) (18) 



is a Boltzmann - like distribution depending on r and t 
only through c(r, t), 

^Eiength = E{r,'L) — Emin 

= XxL'^ + XyL'^ + LLL>c(r),L' = L - L("""). (19) 



where xo = -D^ti^hf^ hf^ and 



i(o) 



i(o) 

V 



y 

Lt^ — Jcm/Xx, 
= — Jcm/Xy 



(22) 



correspond to L 



(min) 
x{y) 



from ([201) provided we neglect 



chemotaxis. 

Eq. (pij) together with Q form a closed set which 
coincides with the classical Keller-Segel system Q. It 
has a finite time singularity (collapse) and was exten- 
sively used for modeling aggregation of bacterial colonies 
[H, [li|. We show below that near singularity contact 
interactions between cells could prevent collapse. 

Eqs. similar to (PT|) and ^ for the case of contact 
interactions (excluded volume constraint) has been ob- 
tained in mean- field approximation 19] . These Eqs. sig- 
nificantly slow down collapse in comparison with (|2ip and 
([3]) . They still have collapsing solutions if initial density 
is not small. (These Eqs. are applicable only for small 
densities.) 



IV. BEYOND MEAN-FIELD APPROXIMATION 
AND REGULARIZATION OF COLLAPSE 

Main purpose of this paper is to derive macroscopic 
Eqs. which do not have collapse (blow up of solutions 
in finite time) for arbitrary initial densities and are in 
good agreement with microscopic stochastic simulations 
for large cellular densities. This requires one to go beyond 
mean-field approximation. 

We conclude from the previous section that random 
changes of cellular lengths result in random walks of cen- 
ters of mass of cells during the time between cell "col- 
lisions" . Significant simplification in comparison with 
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Eq. ^ is that we have now exphcit dependence on spa- 
tial coordinate r but not on L. We use below notation 
{L^^ , Ly"^ ) for average size of a cell neglecting change of 

that size due to chemotaxis. For CPM {L^\l^^) are 
given by (P^ . If we neglect chemotaxis (i.e. set c — 0) 
then during time between collision cell probability den- 
sity p{r, t) is described by linear diffusion Eq. which fol- 
lows from Eq. (PT|) : 

dtp - D2VIP. (23) 

In a similar way, two-cell probability density piji , r2 , t) 
is described by linear diffusion for two independent vari- 
ables ri, r2: 

dtP^D2{Vl^+Vl_)p. (24) 

p{ri , r2, i) represents probability density of a cells 1 and 2 
having centers of mass at ri and r2 at time t, respectively. 
After making change of variables 

r = ri-r2, R=(ri+r2)/2, (25) 

where variable r describes relative motion of cells and 
variable R describes motion of "center of mass" of two 
cells, Eq. ([M)) takes the following form 

dtP^2D2S/lp+^W\p. (26) 

Each collision involving cell 1 or 2 modifies both p(r, t) 
and p{vi,Y2,t). In other words, it effects random walk of 
each colliding cell. 

We describe first the effect of collisions due to excluded 
volume constraint between cells in ID case. Consider 
a pair of neighboring cells which in ID always remain 
neighbors. We assume that at the time of collision two 
colliding cells have the same size. Generally this is not 
true because sizes of cells continuously fluctuate with 
lengthscale SL^ ^ l/(/3^/^Ay^). However, we assume 

as before that these fluctuation are small \5Lx\ ^ L^x^ 
which justifies our approximation. Collision is defined 
as two cells being in direct contact at given moment of 
time and one of them trying to penetrate into another. 
Collison is prevented by excluded volume constraint. In 
the continuous limit e — > each collision takes infinites- 
imally small time. After collision cells move away from 
each other so they are not in direct contact any more. 
Instead of explicitly describing each of these collisions 
we use an assumption that two cells are identical and 
view each collision as exchange of positions of two cells 
[ill . [4^ . From such point of view cells do not collide 
at all but simply pass through each other. They both 
experience random walk as point objects (cells) without 
collisions according to Eqs. and in the domain 
free from other cells (free domain). (The volume of such 
free domain per cell, which has dimension of length in 
ID, is on average (1 — L'^x\{x,t))/p{x,t).) This means 
that we are considering collective diffusion |25i] of cellular 



density instead of trajectories of individual cells. In con- 
trast, self-diffusion describes mean-square displacement 
of an individual cell as a function of time [2^. This ap- 
proach might be important for describing propagation of 
one cell type through space occupied by cells of another 
type. In this paper we consider only collective diffusion. 

While trajectories of cells in free domain are contin- 
uous, positions of cells in physical space change instan- 
taneously at each collision by for the cell colliding 
from the left and by —L^ for the cell colliding from the 
right. The effective rate of cell diffusion is enhanced as 
free space becomes smaller with growth of cellular vol- 
ume fraction. Let's assume that at the initial time t — to 
centers of mass of two cells are separated by an average 
distance l/p{x, t) and that these two cells collide for the 
first time (meaning that previous collisions of each cell 
involved collision with cells other than these two). This 
yields that their centers of mass are separated by dis- 
tance L^'^ at t = to- li moving reference frame is set at 
one of the cells then other will experience random walk 
with doubled diffusion coefficient 2Z?2 (as seen from Eq. 
((26)) ). where D2 is a diffusion coefficient of each cell in the 
stationary reference frame. Relative motion of two cells 
in moving reference frame corresponds to random walk of 
a point- wise cell with diffusion coefficient 2D2. Consider 
continuous random walk. Number of returns to the initial 
position X = L^x^ during any finite time interval after to is 
infinite meaning that number of collisions between given 
pair of cells in physical space is infinite. However, suc- 
cessive collisions effectively cancel each other since they 

change positions of cells by ±Li°\ What really matters 
is whether total number of collision is even or odd. For 
even number of collisions total contribution of collisions is 
zero. While for odd number of collisions contribution to 
the flux of probability for the left cell (at time to) <x. ii"' 

and for the right cell oc — Because total number of 
collisions is inflnite, the probabilities to have even or odd 
number of collisions are equal to 1/2. While the average 
distance between center of mass of cells is \/p{x, t), the 
average distance between boundaries of two neighboring 
cells in physical domain is 

Ax ^ l/p{x,t) - L^^l (27) 

When separation between surfaces of two cells after colli- 
sion reaches Ax from (p7)) we determine that "extended 
collision" between the pair of cells is over. Namely, two 
cells are not closer to each other than to other neighbor- 
ing cells any more. Therefore, probability of them col- 
liding with each other is not higher any more than prob- 
ability of them colliding with other cells. This extended 
collision includes inflnitely many "elementary" collisions 
but its flnal contribution depends only on whether the 
total number of such collision is even or odd. 

To find average time of extended collision we solve dif- 
fusion Eq. in moving frame 

dtp^ = 2D2^lpm (28) 
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with reflecting boundary condition dxPm (^i"'' ,t) = at 
X — L^x'^ and absorbing boundary condition pmiLx^"^ + 
Ax,t) = at X = Li°^ + Aa:;. Reflecting boundary con- 
dition means that cell does not cross point x — L^x^ . In- 
stead of crossing x — l'^x^ cells exchange positions at each 
collision. Absorbing boundary condition corresponds to 
the " escape point" of cell from extended collision. Initial 
condition is Pmix,to) = 6(x ~ L^x'') which is defined by 
initial zero distance between surfaces of two cells. So- 
lution of the Eq. (|28p with these initial and boundary 
conditions results in calculation of the mean first-passage 
time (escape time) Tesc, which is equal to the extended 
collision time in our case. To find T^sc we use a backward 
Fokker-Planck Eq. 

dt^p(x2,t2\xi,ti) = U' {xi)dxiP{x2,t2\xi,ti) 

-DaVl^p{x2,t2\xi,h) (29) 

(see e.g. [131), where p{x2,t2\xi,ti) is the conditional 
probability density for transition from (xi, ti) to {x2,t2)- 
Dq is the diffusion coefficient and U{x) is an arbitrary 
potential. In ID case U{x) = and Dq — 2D2- 

For convenience of the reader we provide here a deriva- 
tion of Tesc which is similar to the one from Ref. [20j . 
Probability that cell remains inside an interval (0, Ax) 
at time t, provided it was located at x at t — to, is given 
by G{x,t) = p{x' ,t\x,to)dx' . Using the stationarity 
of the random walk p{x',t\x,tQ) = p{x' ^Q\x,tQ — t) we 
obtain from (1^51) that 



dtG = -U'{x)dxG + DoV^G. (30) 
Mean escape time Tesdx) for a cell located at a: at t = io 



IS 



Tescix)^- J {t-to)dtGix,t)dt^ j G{x,t)dt. (31) 

Integrating Eq. (I30p over t from to oo and using initial 
normalization G{x,ta) =^ 1 results in 

- U'{x)dxT,sc{x) + DqVIT,sc{x) = -1. (32) 

Reflecting and absorbing boundary conditions for p result 
in similar boundary conditions for Tf,sc{x) : 



0, Tesc{x) 



x=Ll">+Ax 



0, (33) 



which allows us to solve boundary value problem 
(IMl),® exphcitly: 

T,,,{x)^D^^ j e^v[DQ^U{x')\dx' 

X 

x' 

X j exp [- Do iC/(x")]dx". (34) 



Initial condition implies that T^sci 
yields extended collision time in ID: 



{Axf 

iD2 



Te,,(Li°)) and 



(35) 



Consider mean square displacement of position of each 
cell during extended collision in the stationary frame. 
We use notations (Aa;)^ and (Aa::)2 for mean square dis- 
placement of cells 1 and 2, respectively, in the station- 
ary frame. Using (|25|) we obtain that (A.x)^ + (A.x)2 = 
^^f^ + 2(AX)2, where X = R in ID. Because of the 
symmetry between cells 1 and 2 the following relation 

holds {Ax)l = (Aa;)^ = + {AXf. Center of mass 

of two cells experiences random walk with diffusion co- 
efficient D2/2 (see Eq. ^) and (AXf = Tesc 1D2 = 

i^, (A.)f = (Aa:)i = 

Now recall that with probability 1/2 cells exchange po- 
sitions during extended collision meaning that in expres- 
sion for {Ax)i 2 with probability 1/2 term (Aa:)^ should 

be replaced by (2Li°'' + Aa:)^ resulting in the total mean 
square displacement for cell 1 or 2 



(Ax) 



total 



(Aa;)2 (2il"^ + Ax) 



(36) 



From Eqs. (P7)) . ([551) and we obtain effective (nonlin- 
ear) diffusion coefficient for the cell probability density 



D 



l,ef fective 



i'^^)total 



2T, 



esc.l 



i + (Li'")V 



(1 



4°V)^ 



(37) 



In simulations described below we often use not very large 
number of cells N. In case of finite N each given cell can 
collide with A^ — 1 remaining cells which gives extra factor 
(A^— 1) /N and allows one to rewrite the effective diffusion 
coefficient p7p as follows 



1 



(l-iV-i)(Li''))V 



(1 - (1 - Ar-i)Li%)2 



(38) 



The effective diffusion ([38]) results in ID nonlinear diffu- 
sion Eq.: 



dtp = D2V, 



rl + (l-A^-i)(£<"^)V 
(1 - (1 - Ar-i)Li%)2 

-XoVr • [pVrc(r,i)] 



(39) 



Here we added chemotaxis term from Eq. (|2ip which 
is well justified provided Xc 3> Aa; (see also previous 
section). 

Now consider 2D case. First consider cells with disk- 
shaped form instead of cells of rectangular shapes (i.e. 
cellular shapes fluctuating around a disk) with average 

diameter L^"'. Assume, similar to ID case, that is the 
time of the first collision between two given cells. Pre- 
vious collisions of each of these cells involved cells other 
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then these two. We also average over angles relative to 
the center of one of these two cells which implies that cells 
collide at time t = to with equal probability at every an- 
gle. In that approximation there is rotational symmetry 
in the moving frame and all variables depend only on 
the radial variable r = |r| but do not depend on angu- 
lar variables. Change of variables in Fokker-Planck Eq. 
dtp = 2D2^1p results in 



dtp = -2D2dr- + 2D2VIP. 



(40) 



where p = rp. Eq. is equivalent to the ID Fokker- 
Planck Eq. with potential U{r) = —2D2log{r). Back- 
ward Fokker-Planck Eq. ((29|) with x = r, and the same 
U{r) yields in a way similar to Eqs. ([50]) - ([55)1 . an ex- 
tended collision time (mean escape time) 



T 



esc, 2 



Ax{2L<'°^ + Ax) [L^^^y 



8D2 



log 1 



Ax 



(41) 



We assume again that during extended collision time cells 
on average span entire space, i.e. 



(42) 



Using Eqs. (gT]) and (gll) we obtain effective diffu- 

sion coefficient in 2D (here and below the average disk 

(0)n 



diameter is L = ) 



D 



2^ef fective 



1 + ttL'^p 



1 — ttL^p + TrL^plogiiTL^p) 



(43) 



which after modification to include effect of a finite N 
similar to the one used in ID case, results in 



D 



2,ef fective 



1 + {1- N~'^)ttL^P 



.1 - (1 - N-^)ttL^p+{1 - 7V-i)7rL2plog(7rL2p)J 
Eq. (|44p yields the following nonlinear diffusion Eq. for cells fluctuating around disk-shaped form 

1 + (1 - N-'^)ttL^p 



dtp - D^Vr 



1 - (1 - iV-i)7ri2p + (1 - iV-i)7ri2plog(7rL2p) 



XoVr- [pVrc(r,t)]. 



(44) 



(45) 



Notice that both ID effective diffusion coefficient (|38| and 2D effective diffusion coefficient (l44l) depend only on the 
volume fraction (j) {(p — L^^p in ID and (p — irL^p in 2D). Based on that we propose that the effective diffusion 
coefficient for 2D rectangles also depends only on the volume fraction — L^x'^L^^p, which results in the nonlinear 
diffusion Eq. for cells of rectangular shape 

dtp ^D2V 



i + (i-7v-i)4°)4% 



1 - (1 - 7v-i)Lr4> + (1 - A^-i)i^"^4>iog(ii.''4>) 

I 



XoVr • [pVrc(r,i)], 



(46) 



Numerical simulations in the next section confirm that 
Eq. (I46p agrees very well with the MC simulations of 
microscopic dynamics. 



V. NUMERICAL RESULTS AND 
APPLICATION TO VASCULOGENESIS 



One dimensional eell motion 



In macroscopic limit 
placed by 1 in Eqs. 
Eqs. ©and p. 



N 



> 1 factor (1 ^ 
and jSl), 



- TV^i) is re- 
which yields 



If volume fraction — > 1 then nonlinear diffusion in 
Eqs. Il]) and ([2]) diverges. Thus we do not expect any 
blow up of the solution contrary to Eq. (PT|) . Global 
existence of Eqs. (H]) and ^ together with Eq. ^ can 
be studied in a way similar to [43J . 



Figure [3] demonstrates simulation of ID motion of cells 
represented in the form of fluctuating rods. Numerical 
solution of the continuous model was obtained using a 
pseudo-spectral scheme. For both the CPM simulation 
and numerics of the continuous model we used periodic 
boundary conditions, simulation time was tend = 200 
and values of other parameters were chosen as follows: 
Ar = 1, Lt, ^ Lt^ = 3, = Ay = 1.5, e = 0.01 
Jem — 2, l3 = 15, N = 8, fi = 0. Simulation was per- 
formed on the spatial domain < a; < 100 and the initial 
cell density distribution was p{r, 0) = fcoe~*^*^^~^°-'/^°' 



40 60 

X Location 

FIG. 3: Volume fraction (f> = L''^'' p{x,tf.nd) for the 1" 
motion as a function of 2:. Curve (a): Monte Carlo sii 
of the CPM. Curve (b): solution of the Eq. C 
solution of the Eq. (|47|l . Curve (d) solution of the ] 



□ 



with fco determined by normalization N = 8. ^ [~\ 
shows simulations of the CPM (curve a), numerical solu- 
tions of the macroscopic model without excluded volume 
interactions (PT|) (curve b), macroscopic Eq. (|T7|) (curve 
c), and macroscopic Eq. ([39)1 (curve d). Here Eq. 



dtp = iPaVr 



1 



(1 - (1 - iV-i)Li%)2 



-XoVr- [pVrc(r,t)], (47) 



was derived from the equation of state for the ID hard 
rod fluid which allows one to determine collective 
diffusion coefficient from static structure factor and com- 
pressibility (see e.g. [25l. lisj). 

Figure [3] demonstrates that solution of the Eq. ([551) 
is in much better agreement with CPM than both Eqs. 
(1^ and dUl). While the difference between the CPM 
simulation and solutions of the Eq. pT|) and Eq. ([47|) 
is small but clearly exceeds the error in MC simulations. 
Difference between MC simulation and solution of the 
Eq. ([55)1 is within an accuracy of the MC simulations. 

Difference between CPM and Eq. (|T7)) is due to the 
fact that the equation of state for the ID hard rod fluid 
was calculated in Ref. [43 | from grand canonical partition 
function [i^ while diffusion of cells is a nonequilibrium 
phenomenon resulting in the corrections to the equilib- 
rium partition function. Numerous attempts have been 
made to describe dynamics of interacting Brownian par- 
ticles (see [131 and reference there in). Note also that 
difference between CPM and Eq. (pij) results is not so 
dramatic in ID as in 2D because in ID Keller-Segel model 
does not support collapse p^ . 



FIG. 4: p{x,y,t) in 2D as a function of {x,y) for (a) Monte 
Carlo simulation of CPM, (b) Eq. ((2T| and (c) Eq. (|46)l . 



2. Two dimensional cell motion with chemotaxis 

Figures [4] and [5] demonstrate a very good agreement 
between typical CPM simulation and numerical solution 
of the continuous model Eq. pS)) . Both simulations 
were performed on a square domain < x.y < 100 
over simulation time tend = 400. Simulation parame- 
ters' values are as follows Ar = 1, Lx-, = Lx = 4.4, 
^ Xy ^ 1.5, Jem = 2, /3 = 15, At = 0.1, e = 0.01 
and = 15. Chemical field concentration is chosen 
in the form of c{x,y) = 0.2[1 — e ] and 

does not depend on time. Initial cell density is chosen 

in the form of po{x,y) = kQe~^' ^"rm^ ^1 , where ko 
is a constant that normalizes the integral of the cell den- 
sity to N = 15. Numerical solution of the continuous 
model has been obtained using pseudo-spectral scheme 
with 200 X 200 Fourier modes. A large number of CPM 
simulations (600000) have been run on a parallel com- 
puter cluster to guarantee a representative statistical en- 
semble. 

Numerics for Eq. ((2T|) significantly differs from CPM 
simulations indicating that excluded volume interactions 
are important in the chosen range of values of parameters. 



3. Application to Vasculogenesis 

To test the model we studied effect of the chemical 
production rate on the network formation. Our previ- 
ous results [l^ were limited to relatively small chem- 
ical production rate a < 0.2 in Eq. ([3]) because oth- 
erwise chemotaxis resulted in cellular density which was 
too high for applying mean-field approximation. Here we 
use macroscopic Eqs. (j46l) . and compare numerical 
results with the CPM simulations. Figure [6] shows series 
of simulations with different chemical production rates 
a — 0.5, a — 1.5 and a = 3.0. Simulations start with 
initially dilute populations of cells moving on a substrate 
in a chemotactic field, subject to an excluded volume 
constraint. Both CPM and continuous model simulation 
results indicate that stripe patterns are obtained for high 
chemical production rates. Higher chemical production 
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20 40 60 80 100 

Y Location 

FIG. 5: Cross sections of the volume fraction (j> = 
L''°^L^y^p{xo,y,t) in 2D for (a) CPM simulation, (b) Eq. ^ 
and (c) Eq. (1461) . All cross sections are at the same position: 
xo = 57.75. 

rate, by strengthening the chemotaxis and cell aggrega- 
tion process, eventually leads to higher pattern density 
with smaller average distance between two neighboring 
stripes. Structures of the resulting networks obtained 
using discrete and continuous models are very similar to 
each other as well as to the one obtained experimentally 
for a population of endothelial cells cultured on a Ma- 
trigel film 

VI. SUMMARY AND DISCUSSION 

We have derived macroscopic continuous Eqs. (P) and 
coupled with Eq. ^ for describing evolution of cellu- 
lar density in the chemical field, from microscopic cellu- 
lar dynamics. Microscopic cellular model includes many 
individual cells moving on a substrate by means of ran- 
dom fluctuations of their shapes, chemotactic and con- 
tact cell-cell interactions. Contrary to classical Keller- 
Segel model, solutions of the obtained Eqs. ([I]), (O © 
do not collapse in finite time and can be used even when 
relative volume occupied by cells (f> is quite large. This 
makes them much more biologically relevant then earlier 
introduced systems. We compared simulations of macro- 
scopic Eqs. with Monte Carlo simulations of microscopic 
cellular dynamics for the CPM and demonstrated a very 



good agreement for (f> ~ 0.3. For larger density we expect 
transition to glass state [l^. It was demonstrated that 
combination of the CPM and derived continuous model 
can be applied to studying network formation in early 
vasculogenesis. We are currently working on an impor- 
tant problem in vasculogenesis of simulating self- diffusion 
[25t of one type of cells through dense population of other 
types of cells. 

This work was partially supported by NSF grants DMS 
CPM Continuous 

a=0.5 I 




FIG. 6: Simulation of early vascular network formation with 
different chemical production rates a. Ar = 1, Lt^ — Lry = 
0.6, A, = Aj; = 1.5, Jem = 0.002, f3 ^ 15, fi ^ -0.1, Dc = 0.5, 
7 = 0.014, Ate = e^At = 0.01, e = 0.1, te„d = 60. In the 
Monte Carlo CPM simulation A'^ = 15000 cells were randomly 
distributed in a domain < a;, y < 100 with initial chemical 
field at zero. In the continuous model, a uniform initial cell 
density distribution with 5% random fluctuation was used. 
Scale bar is in units of the volume fraction <j!>. 
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